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ABSTRACT 

Searching for dispersed radio pulses in interferometric data is of great scientific in- 
terest, but poses a formidable computational burden. Here we present two efficient, new 
antenna-coherent solutions: The Chirpolator and The Chimageator. We describe the 
equations governing both techniques and propose a number of novel optimisations. We 
compare the implementation costs of our techniques with classical methods using three 
criteria: the operations rates (1) before and (2) after the integrate-and-dump stage, 
and (3) the data rate directly after the integrate-and-dump stage. When compared 
with classical methods, our techniques excel in the regime of sparse arrays, where they 
both require substantially lower data rates, and The Chirpolator requires a much lower 
post-integrator operations rate. In general, our techniques require more pre-integrator 
operations than the classical ones. We argue that the data and operations rates required 
by our techniques are better matched to future supercomputer architectures, where the 
arithmetic capability is outstripping the bandwidth capability. Our techniques are, 
therefore, viable candidates for deploying on future interferometers such as the Square 
Kilometer Array. 

Subject headings: Techniques: interferometric - Methods: observational - Stars: pul- 
sars: general - 

1. Introduction 

1.1. Scientific Motivation 

Studying the high time resolution radio sky has illuminated the physics of our Galaxy, en- 
abled exquisite measurements of physics at the extremes of gravity, density and magnetic field 
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and uncovered a plethora of exotic objects. The main class of object enabling these measure- 
ments is the radio pulsar: a rapidly rotating magnetic neutron star which emits periodic, short 
pulses at radio frequencies. While pulsars are interesting astrophysical laboratories in their own 
right, they can also b e used to test predict ions of General Relativity through observations of single 
pulsar systems (e.g. iKramer et al.l . l2006l ) , search for gravitational waves with groups of pulsars 
(jYardley et al.l. l201Cl and refereri ces therein) , and test theories of matter at the most extreme den- 



sities (e.g. iDemorest et al.l . 120101 ). The short pulses emitted by pulsars undergo propagation effects 
during their passage through the Galactic Inte rstellar Mediurn (ISM ) , which enables measurements 
of the Galactic magn etic field structure (e.g. IVan Eck et al.ll201ll ) and free electron density (e.g. 
Cordes Lazidl2002l) . The re are also many more pulsars to be found: barely 2000 of the estimated 
30000 (jLorimer et al.ll2006l ) have been detected. 



Radio pulsars are chiefly discovered by searching for periodic, dispersed radio emission. By 
contrast, searches for single pulses of radio emission have uncovered other types of objects, of which 



the m ost widely accepted are the so-called Rotating Radio Transients (RRATs) ([McLaughlin et al 



20061 ). RRATs, like pulsars, are rotating neutron stars but emit only sporadicall y, and are being 
studi ed because they may hold the key to the so-called 'missing supernova problem' (jKeane &: Kramer 
20081 ). Searches for single pulses have also yielded a number of intriguing short-duration radio 



transients, which do not f i t the classical mod els of pulsars or RRATs (e.g. iLorimer et al. 



2007 



Burke-Spolaor et al 



2011 



Keane et al. 



2011 



In spite of these discoveries, there is still much to d o, as the the parameter space of radio 
transients is relatively poorly explored (ICordes et al.ll2004l ). Exploring this parameter space opens 
the potential for discovering new object s and physics. Th ese motivations are behind at least eight 



of planning (e.g. iMacquart et al.l . l2010l ). 



ongoing pulsar a nd single-pulse surveys (jMcLaughlinll201ll ). and more surveys are in the late stages 



1.2. Improving pulsar and single-pulse surveys 

When surveying for pulsars and single-pulse sources, a desirable figure of merit is the product 
of instantaneous sensitivity and field of view, known as 'survey speed'. Improving survey speed has 
three important consequences: (1) it reduces the integration time required to reach a flux density 
limit for periodic sources; (2) it reduces the computational requirements to search for pulsars in tight 
binary system^]; (3) it enables a deeper search of the parameter space for single-pulse transients, 
in terms of rarity or faintness. 

In recent years, the optimal approach for maximizing survey speed has been to use large 
steerable and in-earth single dish telescopes feeding multi-beam receivers and wide-band electronics. 
This approach, however, is probably nearing its limit. Steerable single dish engineering has reached 
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its practical limit at diameters of about 100 m and there are limited sites available for large in-earth 
reflectors, which also suffer from a limited view of the sky. Similarly, modern digital electronics 
have improved to the state where 2 GHz bandwidth is readily achievable, but the physics of pulsar 
emission and the interstellar medium limit the majority of the pulsar energy to the range 0.1 
to 10 GHz, so an increase in processing bandwidth is unable to produce large improvements in 
sensitivity. Finally, the field of view of multi-beam receivers cannot grow indefinitely, as a large 
multi-beam receiver simply blocks too much of the dish aperture to be efficient. 

The likely way forward for improvements in survey speed, therefore, is to use arrays of antennas. 
Using a large number of small antennas achieves simultaneously large field of view and sensitivity, 
and therefore survey speed. Using an array does have a dramatic cost, however: for configurations 
of interest to future surveys, the com putational requirem ents increase to the point where the data 



processing becomes almost infeasible (jSmits et al.ll2009l ). 



1.3. The problem: processing requirements 



The desire to use arrays of antennas to simultaneously obtain the large field of view and 
high sensitivity presents major data processing challenges, both in terms of the required data 
and operations rates. These challenges have inspired a number of novel approaches. For example, 
Daishido et al.l (j200d ) proposed the so called Fast Fourier Tr ansform T e lescop e (FFTT) for a pulsar 
survey, which used FFT beamforming, first proposed by IWilliama (jl968l ). and a square array 
geometry. This approach, also known as the Direct Imaging, has attractive properties in terms 
of operations rates and has been extended, with p articular emphasis on 2 1 cm tomography, to 
arrays of regular, ar bitrary h i erarc hies of grids by iTegmark &: Zaldarriaeal (120101) and arbitrary 
array geometries by iMoraled (|2008l ). In a novel experiment, iJanssen et al.l (|2009l ) enhanced the 
field of view over standard techniques, by employing a uniform linear geometry and phased array 
beamforming, w hich introduced d[ eliberate ambiguities into the synthesized beam (i.e. grating 
lobes). Recently, iTrott et al.l (|201ll ) proposed a method for searching for transient sources directly 
visibility space, rather than image-space. The visibility space approach is promising for arrays with 
sparse, arbitrary geometries, as the number of visibilities is much smaller than the number of pixels. 
But, is not yet clear whether this approach will achieve substantial computational savings. 

Some of the above approaches rely on having control over the array geometry. In many cases, 
controlling the geometry is not possible because it is driven by other requirements e.g. the uv- 
coverage. The Square Kilometer Array (SKA) and its pathfinders fall into this category. In such 
cases, one can reduce the processing requirements by falling back to the reduced sensitivity of 
i ncoherent proces sing, analyzing a smaller field of view than available from the primary beam 



(jMacquart 



D'Addario l2010l l or pointing all antennas at different parts of the sky in a "Fly's Eye" mode 



201 ll ). Nonetheless, a fully-coherent, wide-field and computationally tractable system 



is a desirable goal. 
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We also note that the processing requirements are not hmited to the number of arithmetic 
operations. In fact, in modern supercomputing problems, the processing bottlene ck is not the 



numb er of arithmetic operations, but rather the data bandwidth into processor (jLeback et al, 



20081 ). The bandwidth bottleneck has been identified as a key problem for correlators for larg e 



interferometers, with a number of proposed solutions (e.g. iLutomirski et al.l . l201ll . ICarlsonl . |2010| ). 

To our knowledge, the data requirements of classical fast transient detection techniques has not 
been discussed explicitly in the literature, so we consider them in this paper. 



1.4. Two new techniques 

It is in this context that we propose two new techniques, which we have named after the 
frequency swept signals on which they operate ('chirps'). The first technique, which we call 'The 
Chirpolatorl^ operates by correlating the chirps received by pairs of antennas. The second tech- 
nique, which we call 'The ChimageatorH operates by gridding the cross-multiplied voltages from 
all telescopes to form an image at every sampling time. Both techniques are applicable to arbitrary 
array configurations, exploit the full sensitivity of the telescope, and have substantially lower data 
rate requirements than classical coherent techniques. Thus, these new techniques may be favored 
over classical techniques in many regimes of computer economics and array geometry. 

This paper is organized as follows: in section 2 we provide background of the problem of 
searching for dispersion emission in interferometric data, and the classical solutions. In section 3 we 
describe The Chirpolator and section 4 we describe The Chimageator. We describe a simple model 
to compare our techniques with classical results in section 5, and the results of this comparison in 
section 6. We discuss the implications of these two algorithms on future telescope design and science 
outcomes in section 7 and draw our conclusions in section 8. In the appendices we present detailed 
analysis of the algorithms, a number of novel optimizations, and a discussion of implementation 
considerations 



2. Background 

2.1. Dispersion in the Interstellar Medium 

Before being received by telescopes on Earth, electromagnetic waves from an astronomical 
source must pass through the interstellar medium (ISM), a plasma containing non-relativistic un- 
bound electrons. As it travels through the ISM, the wave undergoes dispersion, or frequency- 



portmanteau of 'chirp' and 'correlator'. 

portmanteau of 'c hirp', 'image' and 'correlator'. The 'imaging' part is inspired by the Direct Imaging approach 
of Daishido et al. I (|2OO0l ). 
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dependent delay between two frequencies i^i and V2 according to the following formula: 

t = -^DM [u:,^ - u^^) (1) 
where t is the time from the beginning of the pulse, v\> v^, and the physical constant is: 

II = ~ 4.15ms (2) 

27rmeC 

for frequencies in GHz. The dispersion measure (DM) describes the number of electrons between 
the observer and the emitting source, defined as: 



DM = / riedl (3) 

^0 

where He is the electron density and d is the distance to the source. DM is usually quoted in units 
of cm^^pc. 

If a narrow pulse is emitted by a source, the frequency of the signal received at Earth will be 
exactly the form of the dispersion. For our analysis, we will consider the form of this dispersed 
pulse as a complex voltage time series. To form this time series, we will compute the instantaneous 
phase of the signal, by integrating the instantaneous frequency of the signal. The instantaneous 
frequency can be found by rearranging Equation [H and we approximate it with a Taylor series 
around t = T/2, which yields: 



oo 

= ^a.(t-r/2r (5) 

where V2{'t) is the instantaneous frequency of the signal, t is the time from the beginning of the 
pulse, and T is the time taken for the pulse to traverse the bandwidth of interest B = vi — U2. The 
first three Taylor coefficients are: 



ao = (6) 
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where, 



a 



(9) 
(10) 



A plot of the true dispersion law, and the linear and second order approximations is shown in 
Figure [H 

To obtain the formula for the voltage time series of the signal received on Earth, we first write 
the phase of the signal, as given by: 



27r I V2{t')dt' 
Jo 



1/2 



oo 



2tt 



i=0 

T 



(-00 + ai/2) + t(ao - aiT/A) + fai/2 



(11) 
(12) 

(13) 

(14) 



where we have expanded out the phase to the i = 1 Taylor term. Finally, we can write the complex 
voltage as: 



s{t) = exp(j>(t)). 



(15) 



In a typical telescope system, the absolute phase of the voltage (i.e. the constant term in 
Equation is not important, and the fixed frequency {t term) is removed by down-conversion, 
so only the and higher terms are relevant. Therefore, in the main text we approximate the 
dispersion with a signal for the form: 



s{t) = expiTTjft^) (16) 

where / = ai ~ —B/T is the gradient of the linear frequency trajectory as shown in Figured) 
Signals of this form are known as complex linear chirps. In the appendices we consider the higher 
order terms. 
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928.6 

Delay (ms) 

Fig. 1. — Delay vs frequency for a pulse with DM of 100cm~'^pc, over a frequency range of 
400 MHz centered at 1.4 GHz. The true dispersion (Equation [T]) and the first (linear chirp), second 
and third order Taylor series approximations around the delay midpoint are shown. Third order 
approximation is barely visible as it very closely matches the true dispersion. 
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2.2. A taxonomy of methods 

An astronomer wishing to search for, or study short duration radio pulses, may employ any 
one of a wide range of dedispersion and array processing techniques, as shown in Fig. [2j General 
properties of these methods are shown in Table [TJ We give more detailed descriptions of the classical 
methods in the following sections. 

Table 1: Methods of searching for pulsars and scaling relations for various figures of merit. The 
columns are: the name of the method, the resolution in radians, the number of pixels in an image, 
and the scaling relation for the sensitivity with number of antennas M. Parabolic dishes and equal 
maximum baselines in the u and v directions have been assumed. 



Method 


M 
(rad) 


iVpix 


Sensitivity Scaling 


Power beams 


1.17-^ (primary beam) 


1 


oc Ml/2 


Fourier imaging / Direct beamforming 


\D 

^max 


( ^max \ ^ 

I D J 


oc M 


Chirpolator / Chimageator 


0.844rj;^^ 


;L 92 (fi^'ngxA^a 


oc M 



2.3. Dedispersion 

The effect of dispersion on a short-duration pulse is to smear it out in time. In order to 
determine the emitted pulse shape, or to detect the pulse with maximum signal to noise ratio, 
the effect of the dispersion must be undone, in a process known as dedispersion. Two methods of 
dedispersion can be employed, as described below. 

2.3.1. Incoherent dedispersion 

Incoherent dedispersion involves two steps. First, the raw telescope voltages are passed through 
an analysis filterbank, such as an analog filterbank. Fast Fourier Transform (FFT) or polyphase 
filterbank to form a set of channelized outputs. Each filter output is then squared, and integrated 
over an interval (typically 0.1-lOms) to form a spectrogram, or time-frequency plane. In the second 
step, a range of trial DMs are searched, by summing across frequency channels, after delaying each 
frequency channel according to the trial DM of interest. We call this method 'frequency incoher- 
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ent' because only the filterbank amplitudes are summed, and the phase information is discarded. 
Incoherent dedispersion is typically used in pulsar surveys because the filtering can be computed 
only once and a range of DM trials can be performed on the same filtered output relatively cheaply. 

While early workers used analog filterbanks, more recent projects digitally sample ba seband 



voltag e signals and perform the digital filtering and dedispersion in hardware and software. I Taylor 



(|l974l ) proposed a computationally efficient method of forming performing incoherent dedispersion, 
known as the 'tree' method, which requires fewer additions than a naive implementation, but 
assumes linear dispersion. The linear assumption can be relaxed by adding padding channels, 
which marginally increases the computational cost. 



2.3.2. Coherent dedispersion 

Coherent dedispersion operates on raw telescope voltages, and involves convolving the telescope 
voltages with the impulse response corresponding to the inverse of the ISM (i.e. the inverse of 
Eq. [T]), thereby forming the maximum signal-to-noise ratio filter, or 'matched' filter. We call 
this method 'frequency coherent' as the data are processed without discarding the phase. Coherent 
dedispersion is used during pulsar monitoring, when the DM is approximately known, as the inverse 
filtering preserves the emitted pulse shape more faithfully than incoherent dedispersion. Coherent 
dedispersion is not used for pulsar surveys, as the computational cost of performing multiple DM 
trials is prohibitive. Coherent dedispersion is most often performed on digitally sampled complex 
baseband signals. 



2.4. Array processing 

When using an array, the question arises of how to best combine the signal from two or more 
antennas. In this section we describe three common approaches. 



2.4- 1. Power beams 

Power beams are formed by envelope detecting the output of each antenna, and summing 
the resulting powers across antennas. We call this method 'antenna incoherent' as the envelope 
detection removes the phase information before the sum across antennas. The power beam is 
sensitive to the entire sky, as long as the integration time of the envelope detector is longer than 
the largest geometric delay, and is usually limited by the primary beam of the telescope antennas. 
The penalty for power beams is that the sensitivity is poor, as it scales as M^/^, where M is the 
number of antennas. The output of the power beams can only be incoherently dedispersed, as the 
phase information is discarded by the envelope detector at the antenna. 
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2.^.2. Direct beamforming (Tied array beams) 

Direct beamforming involves delaying the voltage signal from each telescope to compensate 
for the array geometry and summing the resulting voltages. This technique is also known as 'tied 
array beam forming'. This method 'antenna coherent' as the phase information is preserved. The 
resulting beam has the size of a synthesized beam which is much smaller than the telescope primary 
beam. Unlike power beams, the full array sensitivity is preserved as it scales with M. As a tied 
array beam provides a voltage stream, either coherent or incoherent dedispersion can be used. 
Multiple tied array beams can be deployed to increase the field of view. 



2.4-3. Fourier imaging 



Fourier imaging involves cross-correlating the telescope voltages with one another to form a 
set of complex 'visibilities' which are Fourier transformed to form an image. Cross correlation can 
be performed either by an initial filtering step followed by cross-multiplication (so called FX corre- 
lation), or cross correlation followed by a Fourier transform (so called XF correlation). Each pixel 
of the image must be separately incoherently dedispersed, as the pixels are spatially independent. 
Coherent dedispersion cannot be used because each pixel contains only amplitude information. 

Fourier imaging achieves the full array sensitivity over the full primary beam of the individual 
antennas. As only incoherent dedispersion can be used, Fourier imaging is most suited to surveys. 
Fourier imaging requires a so called 'corner-turn', or matrix transpose between the imaging and 
dedispersion stages, which can result i n very hig h data rates betweeii the two steps. Fourier imaging 
has recently been used by iLaw et al.l (120111 ) and IWayth et al.l (|201ll ) to giant pulses from the Crab 
pulsar. 



3. The Chirpolator 

In this section we provide an intuitive description of The Chirpolator, and provide a derivation 
of the equations beginning with the simplest two antenna case. We then extend the results to 
multiple antennas in ID. Extensions to 3D telescope geometries, non-linear dispersion delay and 
novel techniques for efficiently implementing The Chirpolator are described in Appendix [Xl 



3.1. Intuitive Description 

Here we describe an overview of The Chirpolator to aid the intuition of the reader. Put 
simply. The Chirpolator exploits the observation that when a linear chirp received by one antenna, 
is multiplied by a delayed linear chirp received at another antenna, the result is a fixed-frequency 
tone whose frequency is proportional to the geometric delay. The DFT of these tones can be 
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coherently combined across all antenna pairs to form a detection metric. 

A more rigorous mathematical description is described in Section 13.21 and following. 

1. We model the dispersed pulse from an astronomical source as a finite-duration linear chirp 

i.e. a signal whose frequency sweeps linearly across the bandwidth (B) in a time T 
(see Figure [3l top panel). Such a signal has a constant frequency gradient of / = B/T. 

2. This signal is received by two antennas, indexed p and q. The signal at antenna pi s 
delayed with respect to the signal {sq{t)) at antenna q by an unknown geometric delay (r). 

3. The difference in frequency between the two signals is constant for the duration of the pulse, 
and is equal to /r (Figure [3l top panel). 

4. If we multiply the signal from antenna p by the conjugate of the signal from antenna q, the 
result {xpq{t)) is a tone at fixed frequency (Figure El bottom-left panel). This multiplication 
is equivalent to 'downconversion' (also known as 'mixing'), which shifts the frequency of a 
signal in a radio frequency system. In the mixing case, an incoming signal is multiplied by 
a fixed- frequency Local Oscillator (LO), and the result has a center frequency which is the 
difference between the center frequency of the incoming signal, and the LO frequency. In 
our case, both the 'LO' and the incoming signal are sweeping at the same rate (/) but the 
frequency difference remains fixed. Thus, the signal at antenna p is effectively 'downconverted' 
by an 'LO' (provided by antenna q) which is perfectly matched in frequency, yielding a fixed- 
frequency tone. 

5. We have a fixed frequency tone, with unknown frequency (/r) and duration T. In practise, 
this tone will also be contaminated by noise. By taking the Discrete Fourier Transform (DFT) 
of this signal (Xpg[A;]), all the energy of the sinusoid is coherently added into a small number 
of DFT bins, while the noise adds incoherently. Therefore, taking the DFT increases the 
signal-to-noise ratio by approximately the square root of the DFT length. Also, for most 
arrays of interest, the signal can be more compactly expressed in a DFT as the range of 
possible frequencies is much smaller than the number of samples (see Section [A.3.ip . which 
reduces the downstream data and processing rates. 

6. We repeat the above two steps for each antenna pair, and produce a DFT spectrum for each 
(Figure [3l bottom-right panel). The spectrum has a peak at frequency /cq- This frequency is 
proportional to the geometric delay (r), which is in turn proportional to the baseline length 
(bpg) and angle of arrival (6) (see Figure H]). The value of the peak of the spectrum {Xpq[ko]) 
is a complex number whose phase {^pq) is also a function of the geometric delay (r). 

7. Finally we form an image, which is a detection metric {P{9)) for a range of trial directions 
of interest. To produce the detection metric for a given direction of interest, we compute the 
expected arrival frequency (feg) and DFT phase correction {^Iq) for a given antenna pair. We 
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then pick out the DFT bin at the expected frequency {Xpq[k'Q]) and multiply by the phase 
correction {^pg) so that the bins for all pairs have the same absolute phase (See Figure [7|)- A 
vector sum of the phase-corrected DFT bins over all antenna pairs is a coherent sum across 
all antennas, and yields a detection metric in the direction of interest. 

8. In practise, both the time of arrival, and actual DM (equivalent to T) are not known in 
advance. Therefore, we repeat the above procedure in a sliding window fashion, and assuming 
a range of DMs. This repetition can be efficiently implemented using a number novel of 
techniques (see Section [A. 3p . 



3.2. Two-antenna case 

In this section we develop a more rigorous description of The Chirpolator. We begin by 
considering a single pulse which has been dispersed by the ISM, which we approximate by a linear 
chirp impinging on an ideal (perfectly calibrated) two-antenna array. A schematic of the scheme is 
shown in Fig. [3l 



This technique was described by iGershman et al.l (l200ll ) as the maximum likelihood detector 



for a single chirp, which they termed a 'chirp beamformer'. 

As described in Section 12.11 the voltage waveform received by an antenna can be written as a 
complex linear chirp with unit amplitude: 

s{t) = exp (TTjft^) (17) 



where / = —B/T is the chirp rate, B is the system bandwidth and T is the time taken for the chirp 
to cross the bandwidth. Assume this signal is received by two antennas, with the signal delayed at 
antenna g by r seconds with respect to the arrival at antenna p. The product of the chirp received 
by antenna p, with its delayed and conjugated counterpart from antenna q is: 



xp,{t) = sp{t)sl{t) (18) 
= s{t)s*{t-T) (19) 
= exp (vrj/t^) exp (-7rj7(t - r)2) (20) 

= exp (7rj7(2tr - r^)) (21) 

which is a complex sinusoid of frequency /r and phase — vrj/r^. Taking the product in this way is 
also termed 'mixing'. We have assumed here that r <C T, which implies that the signals received by 
both antennas substantially overlap in time. If we have discrete-time sampling we simply replace 
t —7- n/fs where n is the sample number and fs is the sampling frequency. For complex Nyquist 
sampling fs = B. 



The sampled version of Xpq{t) is, therefore 

Xpq[n] = exp (^TTjf - r^^ ^ 



If we take the Discrete Fourier Transform (DFT) of Xpq[n] over N samples, where N = 
BT, and by using the standard result of the DFT of a complex sinusoid of finite duration 
obtain: 



Xpq[k] = DFT{xpg[n]} 

^ f -27r jnk \ ( ,./2nr ^\\ 

( -2'Kjnk 2-KjnfT\ / . ■ 2\ 
= N + f """"P [-^3fT ) 

n=0 \ is J 

= exp (-TTj/r^) exp (^-^^ (^k - ^ 

= ^pq{k-ko)DN{k-ko), 

where fco the frequency of Xpg [k] (in units of DFT bins) , given by 



{B/T)T{fsT) 

fs 

= Bt, 

D]\f(x) is a real- valued amplitude term, whose shape is the Dirichlet kernel, defined as: 

( N x = 

DN{x) = i x^O,-N<x<N 

sm(-ivx/N) I ' 

and $pg(x) is a unit-amplitude complex phase term given by: 



^Use the shift theorem, the Fourier transform of a delta function and the similarity theorem. 
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^pg{x) = exp (-TTjfT^) exp { -irjx — ) (32) 



w 

~ exp {—njx) . (33) 
For the two antenna case, we can write the geometric delay r simply as: 

r = (34) 
c 

where bpg is the distance between the antennas, 6 is the angle of the source off the phase center, 
and c is the speed of light in the medium. The discrete frequency of Xpg, k^, corresponds to the 
position of the peak in the spectrum Xpq[A;], and is related to the baseline length, angle of arrival 
and bandwidth by: 



ko = Bt 

_ ^ bpg sin e 
c 

Thus, the frequency of the mixed signal is linearly related to the sin of the angle arrival and 
the baseline length, as sketched in Figure [H 

To get an idea of the important factors in the above expressions, we can substitute typical 
values for current medium-sized dish-based radio telescopes. For a baseline of 1 km, bandwidth of 
400 MHz, dispersion delay of 1 s and beamwidth of 1 degree, <C 1 and the first exponential term 
of Equation 1331 approaches 1. In the same regime, the number of samples is large, of the order 
N = fsT = BT > 10^. The large number of samples has two consequences. Firstly, in Equation 
[33l (A^ — l)/N — 1, and, most importantly the amplitude term, D]\f{x) in Equation [31] has only a 
small region of support around /cq- This allows computational savings by allowing us to truncate 
the computation of DFT bins to a few bins centered around k^, meaning that calculation of the 
full DFT spectrum is not required, and downstream processing is also considerably reduced. 



(35) 
(36) 



3.3. Multiple telescopes in ID 

To detect a chirp with a given / coming from an unknown direction, we form a detection 
metric, or intensity image, over a range of directions of interest. The detection metric is formed 
by phasing up results from all pairs of antennas. For simplicity we assume the array is perfectly 
calibrated, we ignore the smearing from the D]\i[x) term. Assuming the DFT spectrum is a single 
delta function, with all energy in the bin: 
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k'f) = round(ko). (37) 

then for a particular direction of interest 6, we compute k'^ for each antenna pair using Equations 
[36] and [371 and then compute the value of a single DFT bin Xpg [k'^] . We can essentially phase up 
each antenna pair by multiplying by the conjugate of the known phase term in Equation [33] and 
a detection metric can be formed by performing a vector sum across the phased-up antenna pairs 
according to: 

M-l M-l 
p=0 q=p+l 

The procedure can be repeated for a range of 9. If the direction of interest 9 and the actual 
angle of arrival coincide, the DFT will have a peak at kQ with a value of Xpg[/i;Q]. By substituting 
Equation [27] P{9) reduces to: 



M-l M-l 

^(^) = E E %,{K-ko)<!>M-ko)Dr,{k',-ko) (39) 

p=0 q=p+l 
M-l M-l 

= E E ^^(^o-^o) (40) 

p=0 q=p+l 
M-l M-l 

=^ E E ^^(0) (41) 

P=0 q=p+l 

M(M - 1) 

= ^ (42) 

In general, the quantity P{9) will be complex- valued. We are not interested in the absolute 
phase of the signal, so a more useful metric, for thresholding is: 



m = \P{9)\' (43) 

The fact the sum over antenna pairs is a vector sum means the resulting signal-to-noise scales 
with M rather than the M^/^ scaling for non-coherent addition. 

The time sequence of E{9) for a particular value of 9 can be considered a typical time sequence 
of power measurements, and can be subjected to the usual pulsar detection methods such as 
periodicity and acceleration searches. 
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3.4. Compensating for the smearing in D]\f{x) 

Equations [37] and [39] assume that all the energy is concentrated in a single bin. For an arbitrary 
angle of arrival, this is not the case, and in the worst case, the energy can be spread over all the 
bins in the DFT (see Figure [3]). For large N, D]\f{x) has relatively compact support, so we can 
truncate the number of DFT bins we compute, as well as the number of bins which need to be 
summed for a given direction of arrival and antenna pair. We can choose, therefore, to truncate 
the computations to 2F + 1 bins centered around feg. In practice one can choose a value of F that 
provides the best trade between computational cost and signal-to- noise ratio (SNR). 

To capture the energy with support [—F, F] around ko, we perform a matched filter operation 
against the expected amplitude response function, which is the shifted D]\}{x). Therefore, equation 
[39] can be trivially generalized to: 

Af-l M-l k=+F 

= E E K,(^ + ^o-ko)DN{k-ko)Xpg[k + k'o] (44) 

p=0 q=p+lk=-F 

4. The Chimageator 

We now describe an alternative method for combining signals from multiple telescopes based on 
gridding cross-multiplied voltages. Once again, we begin with an intuitive description and provide 
more mathematical rigour in later sections. 

4.1. Intuitive Description 

Here we describe an overview of The Chimageator to aid the intuition of the reader. The first 
three steps of The Chimageator are exactly the same as The Chirpolator (Section 13. ip . i.e. The 
Chimageator exploits the observation that when a linear chirp received by one antenna is multiplied 
by a delayed linear chirp received at another antenna, the result is a fixed-frequency tone whose 
frequency is proportional to the geometric delay. The difference between the two techniques is how 
the cross-multiplied data are combined: The Chimageator exploits an efficient spatial FFT at each 
sample time. The result is a dynamic spectrum where the chirp deposits energy along a linear 
trajectory in the time-spatial frequency plane. The gradient of the trajectory is proportional to the 
geometric delay. We sum DFT bins along a range of trial trajectories to form a detection metric. 

1. We begin as with The Chirpolator, by assuming a linear chirp which sweeps across the 
bandwidth {B) in time T, with gradient / = B/T. 
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2. As with The Chirpolator, the chirp is received by two antennas and multiphed together 
(mixed). 

3. The resulting mixed signal {xpq{t)) has constant frequency. Once again the frequency is 
proportional to the distance between antennas, and the angle of arrival. 

4. We would like to take a spatial FFT of the mixed signals over all antennas at each sample 
time. Much like the Fourier transform in regular interferometry, this spatial FFT requires the 
signals to be sampled on a regular grid. To form a regular grid (xj[n]), we take the sampled, 
mixed signal from each pair of antennas (a^pqin.]) and average those products which have the 
same inter-antenna spacing (I), and therefore the same (and therefore redundant) geometric 
delays. This process is known as 'gridding'. Gridding can also be used to interpolate a 
non-uniform array geometry onto a uniform grid so that the FFT can be used. 

5. The gridded signals from all antennas are comprise sort of 'space-time tone'. I.e. for a given 
sample number (n), the spatial frequency of the tone is proportional to the angle of arrival 
(6). Similarly, for a given inter-antenna spacing (/), the temporal frequency is proportional 

to e. 

6. For each sample number n, we take the DFT of the gridded signals over the spatial dimension 
(which can be implemented as an FFT). The result is DFT of a single tone (Xfc[n]), which 
has a peak at the bin ko. 

7. Unlike with The Chirpolator, the peak in the DFT (/cq) is not a constant. In fact the peak 
increases linearly with the sample number n and is proportional to the arrival direction 9. As 
a result, a pulse of duration T arriving from a direction 9 will trace out a linear trajectory in 
time where it will cross a number of spatial DFT bins (Figure [5]) during its duration. 

8. The angle of arrival and DM (equivalent to pulse duration T) are unknown. Therefore, at 
each sample time, we assume a set of trial angles and durations, which correspond to a set of 
trial trajectories. To form a detection metric P{9) for each angle and duration, we sum along 
the trial trajectory (applying a phase correction <I>* as we go). 

9. Additional optimizations are possible. For example, the shorter trajectories can be calculated 
as partial sums along longer trajectories with the same gradient, and the spatial FFTs can be 
averaged before performing the trajectory sums. These optimizations are discussed in Section 

El 

4.2. Formulation for a uniform linear array 

In this section, we develop a more rigourous description of The Chimageator. 
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Consider a linear, perfectly calibrated array of M antennas, uniformly spaced with inter- 
element spacing L. If a linear chirp impinges on the array, the product of the signals from two 
antennas, indexed p, and q is, therefore, given by: 



Xpqit) = Sp{t)sl{t) (45) 

= s{t — pt)s* {jp — qr) (46) 

= exp [tij f{t - prf^ exp (^-Trjf{t - qrf^ (47) 

= exp(-7rj7(2tr(p-g))+r2(p2_g2)^ (4g) 

~ exp(-2TTjftT{p-q)) (49) 



where 

r = - sin(6') (50) 
c 

As before, we can form the sampled signal by replacing t — )• n/ fg. Next, we combine the values 
of for all baselines with the same spacing I, and the sample number n in a process known as 

gridding. The uniform linear array has redundant spacings which can be combined and weighted 
according to: 



Xi[nl 



M-l-l 
p=0 



(51) 



where I runs from 1 to M — 1 (the auto-correlations are ignored). Xi[n] corresponds to the visibility 
measured by combining all baselines with spacing (/ + 1)D and Wi are weights. Wi = 1 corresponds 
to 'uniform' weighting yielding the maximum signal-to-noise ratio but reduced resolution. Wi = 
1/{N — I) corresponds to natural weighting, yields maximal resolution but reduced signal-to- noise 
ratio. Through suitable choice of weights, an arbitrary array be interpolated onto a regular grid as 
required for the spa tial FFT, by griddi ng with a spatially varying set of weights. The interested 
reader is referred to lTavlor et al.l (jl999l . chapter 7, section 3), for details. 



For a single chirp, we can substitute Equation 09] and assuming natural weighting and a uniform 
linear array, the gridded voltages simplify to: 
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AI-l-l 

x'lin] = Yl Wi exp (^-27r j fnT{p-{p + I)) /fs) (52) 

p=0 

M-l-l 

= exp (27rjfnTl/ fs^ ^ Wi (53) 

p=0 



exp(27TjfnTl/fs) (54) 

(55) 



A spatial discrete Fourier transform of the gridded voltages yields: 



Xk[n] = DFT{x[[n]} (56) 
5:'exp(^)exp(2.,/r/-^) (57) 



A/-1 



(58) 

\ / V I ^ / 

1=0 

= DM{k - ko)<^>{k - ko) (59) 

where DM{f) is the Dirichlet kernel defined earlier, and <I> is a unit-amplitude phase term. This 
spatial DFT can be efficiently implemented as a Fast Fourier transform. 

A chirp crossing a bandwidth B = fg in time T, arriving from angle 6 signal puts power in the 
DFT bin given by: 



71 

ko{n,T,9) = JtM- (60) 

Is 

Ti 

= ifTM (61) 

nML , , 

= -^sin(Q) (62) 

^ ^^sin(0) (63) 



which is very similar to the expression for The Chirpolator described previously, with the key 
difference that in this case, the k<^ term now depends linearly with sample number n rather than 
baseline length h^q. Thus, a chirp signal will appear as power along a diagonal trajectory in n-k 
space, as shown in Figure [5l 
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For a chirp beginning at n = 0, the trajectory ends at the DFT bin given by: 



kend = ko{Tf„T,e) (64) 

= /,^sin(^) (65) 
c 

To form an image, we can sum across the diagonal trajectory in DFT bins and time, applying 
the inverse of the phase term to produce an intensity image for a given dispersion delay, by: 



Pt{G) = Y.'^*(^o{n,T,9))Xkikoin,T,e)) (66) 

n=0 

with the scalar energy computed as in Equation 



5. Method of comparison 

We have described two new antenna-coherent techniques for detecting dispersed pulses with 
interferometers. In this section, we describe our method of comparing our techniques to two existing 
classical techniques with roughly equivalent sensitivity: Fourier imaging, and direct beamforming 
with frequency incoherent processing (see Figure [5]). These classical techniques are described in 
Sections [2X2l [2X31 and [2XT] respectively. 

Ideally, we would like to compare the techniques in terms of the true implementation costs. 
But, evaluating the true implementation cost is complicated by a number of considerations: 

• The choice of survey parameters (e.g. minimum & maximum DM, center frequency). 

• The telescope parameters (e.g. number of antennas, system bandwidth, baseline distribution). 

• The economics of available technologies. 

• The details of the implementation on a given technology. For example, how an algorithm is 
parallelized over a number of processors. 

• The techniques do not yield equivalent sensitivities in certain situations (e.g. Section [7.1.2p . 

• The parametrization of the algorithms themselves. 

To explain the final point further: an implementation of an algorithm requires a set of param- 
eters that affects both the cost and the sensitivity of that implementation (e.g. number of channels 
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for interferometric imaging, or F for The Chirpolator). For each algorithm, the relationship be- 
tween the parameters and sensitivity is complicated, and there is no straightforward way to choose 
realizations that yield equivalent sensitivities for all techniques so that their costs can be compared 
fairly. 



5.1. A simple model for evaluating algorithm cost 

To help illustrate, in very approximate terms, the differences in operations and data rates 
required by different methods, we propose a simple model. In this model, we split each algorithm 
up into two basic functional blocks: the processing required before an integrate-and-dump step, 
and the processing required after it. We also consider the data rate required between the two 
blocks, i.e. immediately after the integrate-and-dump step. We acknowledge that this model does 
not consider very important details of how data is transported within each block, and acknowledge 
that the bandwidth bottlenecks may indeed be within each block, rather than between the two. 
But, the bandwidth requirements inside each block are a strong function of the way the the pro- 
cessing is parallelised inside each block, and quantifying the many different methods for doing this 
parallelization are outside the scope of this paper. 

This functional breakdown applies to the techniques as follows: 



The Chirpolator : The pre-integrator step is the sliding-DFT (Section IA.3.2p . The integrated 
output is a set of DFT results per DM trial. The post integrator step is the imaging per DM 
trial. Detail of the data and operations rates are described in Appendix lAl 

The Chimageator : The pre-integrator steps include gridding and integration to the shortest 
sampling interval. The integrated output is a sequence of partially-averaged images. The 
post-integrator steps include the remaining integration for the full range of DM trials, and 
the imaging. Detail of the data and operations rates are described in Appendix iBl 

Fourier Imaging : The pre-integrator steps include cross-correlation and integration. The inte- 
grated output is the visi bilities. The post-integrator steps includes gridding, FFT and tree 
incoherent dedispersion (jTavlorl Il974l ). For the the bandwidth requirement, we sum both 
the requirements for both the vi sibilities , and the 'corner turn' required for dedispersion. 
Operations rates are described by ICordesI (j 19971 ) 



Direct Beamforming : We form as many tied array beams as required to cover the entire the 
primary beam. The pre-integrator steps include the beam forming and integration. The 
integrated output is a power spectrum per bea m. The post in tegrator step is tree incoherent 



dedispersion. Operations rates are described by ICordesI (|1997l ) 
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5.2. Array, survey and algorithm parameters 



To arrive at concrete values of bandwidth and operations rate, we must define a full set of 
parameters for an array, survey an d each a l gorith m. To motivate our example, we choose a set of 
parameters based on the SKA from ICorded (j 19971 ) , as shown in Table [2j Clearly, the evaluating the 
performance of all techniques as a function of all parameters results in a highly-multidimensional 
dataset. For the sake of simplicity, we leave only one free parameter: the number of antennas in 
the array (M). We let M go from 2 antennas up to 2000 antennas, which covers the range of values 
for SKA and its pathfinders. 



6. Results 

For all but the largest arrays, The Chimageator and The Chirpolator have substantially su- 
perior bandwidth requirements than the classical techniques (Figure [6|). The lower bandwidth 
requirements are achieved because of a difference in timescale that needs to be sampled by the 
integrate-and-dump step: our techniques sample the shortest dispersion delay, while the classical 
techniques sample the shortest dedispersed pulse duration. As a dedispersed pulse can be substan- 
tially shorter than the dispersion delay, the classical techniques must dump their integrators at a 
much higher rate, therefore requiring larger bandwidth between the functional blocks. One addi- 
tional factor worsens the bandwidth requirements for Fourier imaging in particular: below about 
~ 100 antennas the bandwidth dominated by the dedispersion 'corner turn'. 

In terms of post-integrator operations rate. The Chirpolator betters all other techniques up to 
~ 200 antennas. This low rate for small arrays is consequence of both the low input bandwidth, 
and the fact the imaging operates on a per-baseline basis. Above ~ 200 antennas, the Direct 
beamforming method is the clear winner, as the dedispersion cost is fixed by the longest baseline, 
rather than the number of antennas. 

In terms of pre-integrator operations rate, Fourier imaging is clearly the most efficient for all 
array sizes of interest, with our techniques requiring between 2 and 4 orders of magnitude more 
operations for equivalent array sizes. 



7. Discussion 

For any array size, there is no clearly superior algorithm in all measures. The Chirpolator has 
high pre-integrator operations rate, but has good post-integrator and data requirements for small to 
medium arrays. The Chimageator has consistently high a post-integrator operations rate. Fourier 
imaging is computationally attractive but has very high data rate requirements, either due to the 
corner turn in small arrays, or the visibility data rate in large arrays. Direct beamforming has very 
high operations and data rate requirements for small arrays but becomes somewhat competitive 
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Table 2: Parameters used in our example model. The parameters for The Chu'polator and The 
Chimageator are defined in Appendices (Aj and iBl respectively. 



Parameter Value 



Array parameters 




System Bandwidth (MHz) 


400 


Antenna size (m) 


12 


Maximum Baseline (m) 


1000 


Center Frequency (GHz) 


1.4 


Number of polarizations 


2 


Survey parameters 




Minimum DM (cm'^pc) 


10 


Maximum DM (cm'^pc) 


1000 


Fourier imaging & direct beamforming 


Number of frequency channels 


1000 


Number of DM trials 


1000 


Integration time (seconds ) 


10-* 


Bytes per visibility (post correlator) 


2 


Bytes per image pixel 


1 


Chirpolator specific 


DM step (e) 


0.1 


Smearing support size (F) 


1 


Time oversampling (Kt) 


4 


Spatial oversampling (ks) 


1 


Bytes per DFT bin 


2 


Chimageator specific 


DM step (e) 


0.1 


Smearing support size (F) 


1 


Time oversampling (fi:t,o) 


4 


Spatial oversampling («s) 


1 


Operations per grid point 


50 


Bytes per FFT bin 


2 
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for larger ones. The preferred algorithm, therefore, will depend on the details of the array, survey 
and algorithm parameters, and the economics of available computing technologies. 

The economics of computing technology are changing rapidly. The increase in arithmetic 
capability of processors has been well described by Moore's law; that is, the number of transistors 
(and by inference, arithmetic capability) on a chip doubles every 18 months. While this prodigious 
improvement is very welcome for the arithmetic part of the problem, it does not hold for data 
rate, which has traditionally grown much more slowly. We propose that, because the arithmetic 
capability of processors is outstripping the bandwidth capability, our techniques with their superior 
data rate performance, will become more and more favorable as technology progresses in spite 
of their requirements for higher operations rate. Therefore, in the time scale of the SKA, its 
pathfinders, our techniques may be preferred over the classical ones. 



7.1. Further work 

7.1.1. Effect of calibration errors 

In our analysis we have assumed an ideal, perfectly calibrated array, in which all the antenna 
gains are equal and have zero relative phase. In practice, each antenna will have uncalibrated 
errors in gain and phase which will affect the performance of our algorithms. While a detailed 
discussion of the effect of calibration errors is outside the scope of this paper, we present here a 
simple proof that phase errors (which we model as delay errors) in The Chirpolator case, will result 
in decoherence across the array and reduced SNR. 

If we assume the uncalibrated delay error between two antennas is Te„ then we can make the 
substitution r — )• r + Tgrr into Equation [35] to obtain the frequency of the tone after mixing: 



k = S(r + re„) (67) 

= ko + fcerr (68) 

Therefore, a delay error changes the frequency of the mixed signal, and shifts the entire DFT 
spectrum from Xpg [k] to Xpq [k + kerr] ■ The shift in the DFTs reduces the amplitude of the detection 
metric, which is formed by a vector sum of the phase-corrected DFT bins from each antenna pair. 
The detection metric has a maximum value when all the phase-corrected DFT bins have the same 
absolute phase. If an antenna pair contains a delay error, the each phase-corrected DFT bin will 
not have the same absolute phase, and the vector sum will not be over a straight line (Figure [7]), 
resulting in reduced amplitude of the sum. This process can be quantified by substituting ko into 
Equation [38l 
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M-l M~l 

PeAO) = E E K<}i''o-k)X,,[ko] (69) 

p=0 q=p+l 
A/-1 A/-1 

= ^ Yl %q{ko-k)^pq{ko + hrr-ko)DN{k'o + hrr-ko) (70) 

p=0 q=p+l 
M-l M-1 

= ^ ^ exp (j7r(/co - ^)) exp (-j7r(/Co + fcerr - A:o) -DAr(-Brerr) (71) 

p=0 q=p+l 
M-l M-1 

= ^ ^ exp {-jirB Terr) Dn{B Terr) (72) 
p=0 q=p+l 

< ^ideal(^) (73) 

The inequality in Equation [73] is a result of the triangle inequality for vector addition (See 
Figured]), and the fact that Dj\[{x) < Dj\f{0). 

It is clear from this argument that delay errors will result in a reduced detection metric, 
resulting in a drop in signal-to- noise ratio. We leave a quantitative analysis of this effect, and other 
calibration effects for future work. 



7.1.2. Extension to millisecond pulsars 

Both our methods have assumed that a chirp is received in isolation, meaning that during 
the duration T of a chirp, no other chirps are received. This condition is violated for millisecond 
pulsars, which have short periods and can have large DMs. The combination of short period and 
large DM means that a chirp will not have finished traversing the system bandwidth B before a 
subsequent chirp is received. 

We can write the isolated chirp condition for a pulsar with period P as: 



P > T (74) 
> ^iDM{v:^;^ - v:^^). (75) 

If the isolated chirp condition is not satisfied, there are multiple chirps occupying the bandwidth 
at any one time. These additional chirps produce additional mixing products at the multiplication 
steps (i.e. in Equation [18] and Il5]l which appear at frequencies that are outside the frequencies 
searched in the isolated chirp case. If only isolated chirp processing is performed, the energy in the 
additional mixing products is effectively lost, with a resulting loss in SNR. Our techniques will still 
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operate effectively, but the SNR achieved will not be as high as with processed by other methods. 
Quantifying loss of energy to mixing products, and resulting loss in SNR, is outside the scope of 
this paper. 

To determine what fraction of pulsars violate this condition, we use the ATNF pulsar catalog 



( Manchester et al.l l2005l rl. This catalog contains the DM and period for all known pulsars. At a 



bandwidth of 400 MHz at 1.4 GHz, 30% of the known pulsars have periods that are too high to 
satisfy the isolated chirp condition for their DM (Fig. [8]), indicating this effect is important. 



8. Summary 

We have described two new techniques for detecting dispersed pulses with radio interferometers, 
which we call The Chirpolator and The Chimageator. These techniques have antenna-coherent 
sensitivities in the isolated chirp case and substantially lower data rate requirements than other 
coherent methods for realistic array configurations including the SKA and its pathfinders. For small 
to medium array sizes The Chirpolator is also more efficient than classical techniques in terms of 
post-integrator operations rate. While the pre-integrator operations rates our methods high in 
some cases, the changing economies of computer design may flavor lower bandwidth requirements 
of our new techniques in spite of their high operations rate requirements. 

KB acknowledges the support of an Australian Postgraduate Award and a CSIRO top-up 
scholarship. The Centre for All-sky Astrophysics is an Australian Research Council Centre of 
Excellence, funded by grant CE11E0090. We thank the anonymous referee for their prompt and 
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A. The Chirpolator: Analysis and implementation 

In this appendix we describe additional extensions to The Chirpolator to include multiple tele- 
scopes in 3D and non-linear dispersion. We describe novel methods for efficiently implementing The 
Chirpolator and also derive equations for the resolution and data and operations rate requirements. 



A.l. Multiple telescopes in 3D 



The generalization to arbi trary arrays of elem ents in three dimensions is most easily done in 
the notation of interferometry (jTavlor et al.lll999l . Chapter 2). 

If we measure the [u,f,u)]"^ baseline vector in units of distance (not wavelength), then the 
geometric delay for a 3 dimensional array is: 



ul + vm + w ^\/l — l"^ — m"^ — 1^ 



(Al) 



This preprint was prepared with the AAS lATp^X macros v5.2. 
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where I and m are the direction cosines in the u and v directions respectively. / and m define the 
angle of interest analogous to 9 in the ID case. 

The method of computing the intensity image then proceeds in much the same manner, with 
ko computed with equations [35] and lAH and with P{9) evaluated over two angular dimensions 
instead of one. 



In the main text, beginning at Equation 1171 we have assumed a linear chirp. In fact, at most 
frequencies and bandwidths of interest (i.e. below 10 GHz, and bandwidths >100 MHz), the cold 
plasma dispersion law is much more accurately modeled as cx as shown in Equation [1] and 
Figure [TJ In this section we describe the effect of the true dispersion law on Chirpolator processing 
(decoherence) , and propose a solution (oversampling) . 



To determine the effect of the higher order terms on Chirpolator processing, we begin by 
considering the frequency of the mixed signal Xpq (Equation [T8|) , which is the difference between 
the instantaneous frequencies of the signals from the two antennas. Assuming a delay t <^T, the 
instantaneous frequency difference between the two chirps is given by: 



where we have used the Taylor expansion described in Equation [5l From Equation IA3I we can see 
that the effect of nonlinear dispersion on The Chirpolator processing is to smear out the signal across 
a wider range of frequencies after mixing the two antenna signals (Fig. [9]). The departure of the 
frequency from the linear assumption is significant for typical array configurations and dispersion 
(Fig. [TOl) . and is worst far from the phase centre, on the long baselines, and at t = 0, where it 
can be approximated as the difference between the linear approximation and the 3rd order Taylor 
series (Equation IA3P : 



A. 2. Non-linear dispersion 



A. 2.1. The Problem: Decoherence in the DFT bins 





(A2) 
(A3) 



aiT - 2a2T{t - T/2) 




(A4) 



When a signal with non-constant frequency is passed through a DFT, the amplitude of the 
DFT output is reduced, which we call decoherence. 
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We identify three regimes in which the system operates: 

• The smearing is ^ 1 bin, in which case the decoherence is small and can be ignored. 

• The smearing is ~ 1 bin, in which case the signal still occupies only one bin, but the deco- 
herence within that bin is significant. In this case, the DFT must be broken into a number 
of sub-integrations, with each sub-integration requiring a complex phase rotation to recover 
the coherence. 

• The smearing is > 1 bin, in which case there is energy in multiple bins. The DFTs must 
be broken into a number of sub-integrations. The final output must before formed with is a 
complex phase rotation of a range of sub-integrations of different DFT bins. 

If the smearing is > 1 bin, (e.g. Fig. [9]), the true dispersion occupies a higher DFT bin 
than the linear assumption for approximately half the pulse duration. To capture energy from the 
higher frequencies, additional DFTs must be computed that would not be required under the linear 
assumption. In the nonlinear case, the maximum number of DFT bins increcLSGS from ho max 

to 

:tnax ~\~ '^mixi which increases the operations and data rate requirements for the DFT step. At the 
worst case longest baseline of 1 km, at 1.4 GHz, 400 MHz and 0.5 degrees from the phase centre, 
^mix — 7, and the number of DFT bins required increases by a factor of 5mix/fco,max — 

63%. The 

additional DFT bins increases non-linearly as a function of baseline, so accurately estimating the 
total increase over the whole array requires a knowledge of the exact baseline distribution. To 
obtain an approximate figure, assuming a baseline distribution where the mean baseline length is 
half the maximum baseline length, we propose that the total increase is approximately half the 
worst case figure, i.e. 32%. 



The key to handling the nonlinear dispersion, therefore, is to dump the integrator more often 
than required for the nonlinear case (oversample) , and phase-correct the results to obtain the 
coherence again. To quantify the amount of oversampling required where the smearing is > 1 bin 
we need to quantify the response of the DFT to the mixed, non-linearly dispersed signal. As shown 
in FigO the frequency of the mixed signal is well approximated by the third order Taylor expansion 
described in Equation [5l The phase of the mixed signal is, therefore, given by the integral: 



A. 2. 2. The solution: oversampling 




(A5) 



27r{—aiTt — a2Tt{t — T)) 
27r(t(— aiT + a2TT) — t^a2T). 



(A6) 
(A7) 
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The sampled, mixed signal can then be expressed as: 



m[n] = exp(j(/)mix[n]) 
= exp 



2vrj f Y^-aiT + a2TT) - ) 



(A8) 
(A9) 



and we take the DFT over = fgT samples to obtain: 



N-l 



Xm[k] — exp {—2TTjkn/N) m[n] 

n=0 
N-l 

= exp 2TTj [n{-aiT + oztT - k/N) + n'^{-a2T/ f^)] 



(AlO) 
(All) 



n=0 



The term that is hnear with n has already been dealt with in Equation [271 ^iid is simply the 
DFT of a single tone, so we turn our attention to the term and define the sum: 



L-l 



G{a, L) = exp (lirjan^^ 



(A12) 



n=0 



where: 



a = -a2T/f^ 



(A13) 



From Equation IA12I it is clear that G(0, L) = L, and that for non-zero values of a, the av? 
term introduces oscillations, effectively moving the instantaneous frequency into the adjacent DFT 
bins, so that |G(a, L)| < L for non-zero a. 

We want to determine how large L can be made before some fraction of the energy will be 
lost to adjacent DFT bins. Equation IA12I defines the result of summing a chirp with an initial 
instantaneous frequency of zero, which is essentially the centre of the DFT bin. If we define the 
the coherence loss, or loss in amplitude as: 



7^ = \G{a,L)\/L 



(A14) 



Therefore, the value of L that maintains a required r] is the number of samples to traverse 
half the DFT bin and maintain a given loss. To calculate the required oversampling for chirp that 
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crosses and entire DFT bin, we can pose the question: what oversamphng factor Kt = N/2L is 
required to maintain ry(a, L) above a specified threshold? 

Do get an approximation of the required oversampling factor, we have simulated a typical 
case for the SKA, with a case with a one-sided frequency smearing of the order of 7 DFT bins, 
which reasonably large in the context of Figure [TOl We conclude that a 4 times oversampling yields 
r/ ~ 99% (Fig. [TT]) . We have found empirically that the required oversampling is independent of 
DM. 

Because the time of arrival is not known, a one would typically require ~ 4 times oversampling 
to obtain a sample which is integrated over a large fraction of the incoming signal. The equivalence 
of the oversampling rates required for time oversampling, and nonlinear dispersion correction, 
implies that nonlinear dispersion does not substantially drive the oversampling in this instance. 

A. 3. Implementation optimizations 

In the main text we assumed a single value of / (equivalently a single value of the DM) , and that 
the DFT window is exactly time-aligned with the chirp. In practice, neither the time of arrival for 
the chirp, nor the / are known in advance and we would like to maximize our chances of finding the 
signal. The maximum likelihood approach to the problem of maximizing the detection probability 
when the waveform parameters are unknown, is to pass the signal through many different matched 
filters, each with a particular realization of the unknown parameters. In our case, we would evaluate 
P{9) and ^pg[fc] independently on a range of values of / and on a set of overlapping windows in 
time. 

Significant computational savings can be made as described in the following sections. 

A.3.1. Compute only DFTs required on a baseline basis 

We do not have to compute the same number of DFT bins for each pair of antennas. In fact, 
for a given pair of antennas, we only have to compute the DFT for values of k up to approximately 
^o,max) as illustrated in Fig. [H With values from typical radio telescopes, ko^max — 100 for the 
longest baselines, and ko^max — 5 for the shortest baselines. If the baseline distribution is such 
that the average baseline is half the maximum baseline, this strategy saves a factor of 2 in DFT 
operations and data rate. 
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A. 3. 2. Efficient calculation of Xpq[k] with sliding DFTs 

We consider problem of computing DFT values for overlapping time windows. For typical 
array configurations and DMs, the number of samples in the DFT (N) is of order 10^, whereas the 
number of usable DFTs is of the order /co,max — lO'^, meaning that computing a full FFT would 
result in a very large number of unused bins. In addition, we do not require a DFT result every 
sample, which means a sliding window DFT result every L < N samples is adequate. 

A naive method to computing the sliding window DFT is to (1) compute the dot product of N 
input samples with a complex sinusoid of appropriate frequency, then (2) shift the input sequence 
by L < samples, and (3) compute the dot product on the shifted samples, with the same complex 
sinusoid. This naive method requires A^ complex multiplications per L samples, per DFT bin, and 
corresponds to an operations rate of roughly fgN/L per DFT bin. 

Jacobsen &: Lyons (j2003l ) describe a 'the sliding DFT', a more efficient method for computing 



a small number of DFT bins in a sliding window manner. The sliding DFT is a recursive filter that 
produces a sliding window DFT output according to: 



Sk[n] = Sk[n - 1] exp {-2piTTjk/N) - x[n - N] + x[N] (A15) 

where Sk[n] is the sliding window DFT output for sample n and bin k, and x[n] is the sampled 
input sequence. Equation IA15I is effectively a moving average filter implemented as a Cascaded 
Integrator Comb (CIC), with a complex resonator embedded in the integrator feedback path. The 
sliding DFT has a operations rate of only ~ 3/^ per DFT bin, which is significantly less than that 
required for the naive method. 

In practice, we do not require an output every sample, so the operations rate can be further 
reduced by computing a block-based sliding DFT. In this case, we compute the partial DFTs, 
time-indexed by m in blocks of L samples: 

L-l 

^fei*^] = exp (— 27rj(n -|- mL)k/N) x[n + mL] (-^16) 
n=0 

and the and form the the DFT over the full number of samples A^ by applying a moving average 
filter on the partial DFTs: 



Sk[m] = Sk[m-l] + Vk[m\ 
This method is illustrated in Fig. [T2j 



-Vk[m-N/L]. 



(A17) 
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The block-based sliding DFT has a lower operations rate than the sliding DFT, because the 
moving average (CIC) stage (Equation IA17P operates at the block rate, rather than the sample 
rate, which resulting in an operations rate of ~ + 2fs/L = fs(l + 2/L) ~ fs- 



A. 3. 3. Factorizing the DFTs 

As the DM and therefore the value of / is unknown, a search in / is required to maximize 
signal-to-noise ratio of pulse at unknown DM. This search through / is equivalent to varying the 
size of the DFT: A'^. One might choose to use a bank of DFTs, each with a length of A^ = NqcI, 
where Nq is the length of the DFT corresponding to the shortest DM of interest, and d is a positive 
integer. If we have a bank of DFTs, each starting at the most recent sample and extending back 
in time by N samples, then we can factorize some of the DFTs by noting that some of the basis 
functions for the long DFTs can be formed by concatenating the basis functions for the short DFTs. 
By way of example, the result of the 5*2 bin for the length 2Ao window can be trivially computed 
by summing the adjacent, non-overlapping results of 5i over the length Nq window, as illustrated 
in Fig. \M 

If we write 5'^[n] as the DFT result for bin k at sample tim^ n for a length N DFT, we can 
say that a DFT of length Aqc? can be computed from the sum of D shorter length NqcI/D DFTs, 
if it can be written as: 



n 



(A18) 



d'=0 



where k and d are integers. 

The bin 5f "'^ can be factorized if and only if d/D and k/D are integers, that is d and k must 
have a common, non-unity factor D which implies that d and k canno t be co-prime. The pro bability 
of two random integers being co-prime is approximately 61 per cent (jHardv &: Wrightll2008l ). which 
implies that approximately 39 percent of DFT bins can be factorized. If the shorter DFT results 
are already available, computing the factorized DFT requires D — 1 operations, which is trivial 
in comparison with A operations to compute the full DFT. As a result, DFT factorization saves 
approximately 39 percent in complex operations. 



In practice, one would compute the factorized DFTs on the sliding DFT block outputs indexed by m. We have 
kept the full sample rate n here for clarity. 
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A. 3. 4- Efficiently computing negative DFT bins 

Computing the DFT output requires the multiphcation of the complex input sample with 
the complex exponential. On a standard computer, the complex numbers are stored as real and 
imaginary parts, and the complex multiplication is performed in the following way: 

rfc[n] = x[n] exp(-27rjA;n/iV) (A19) 
= {a + jh){c-3d) (A20) 
= ac + hd + j[-ad + bc) (A21) 

where rf[k] is the result of the multiplication of the input sample with the sinusoid of frequency 
/c, a and h are the real and imaginary parts of the complex input sample, and c and d are the real 
and imaginary parts of the complex sinusoid. 

Assuming the phase center is set to the center of the primary beam, the DFTs must be 
computed for frequencies over the range [— /^o.max, ^o,max] to cover the full field of view. Therefore, 
each positive bin has a negative counterpart. To compute the negative frequency, we could also 
separately calculate: 

r_k[n] = x[n]exp{2TTjkn/N) (A22) 
= ia + jb){c + jd) (A23) 
= ac-bd + j{ad + be) ( A24) 

The calculation of rk[n\ and r_fc[n] naively requires 12 operations (8 multiplications and 4 
additions). But the multiplications are common between the two results (Equations IA21I and 
IA24p . which means both results can be computed with only 8 operations (4 multiplications and 4 
additions). This results in a saving of 33 percent over the naive implementation. 

A. 4. Performance 

A. 4-1- Resolution 

From the definitions of Equations [38l and HHl the amplitude of the response to a chirp on a 
given baseline is |P(0)p = D'j^{x). Therefore we approximate the spatial resolution of the The 
Chirpolator as the Full Width Half Maximum (FWHM) of D%{x) on the longest baseline. The 
FWHM of D%{x) is defined by: 
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D%i2xFWRM) = Id\0) (A25) 



sm(27rxFWHM/A^) V2 

Taking the 3rd order Taylor expansion of the sin terms, and solving for the non-trivial solutions 
of X, we obtain: 




1 /6(l-l/V2) 

(A28) 
(A29) 

To convert Equation IA291 which is the width of the main lobe in units of DFT bins, to an 
angle, we rearrange Equation 1361 which yields: 



sin^ = (A30) 

and by applying the small angle formula, and substituting Equation IA29I as the DFT bin resolution 
(i.e. Ako = xfwhm), we obtain the spatial resolution of The Chirpolator: 

AO = Ako-^ (A31) 

= 0.844 (A32) 

-t> Orn a X 



A.4-2. Operations rates 
We compute the number of operations required to form images o f the ful l field of view of a 



telescope comprised of parabolic dishes. We keep to the convention of ICordesI (|l997l ) of counting 
complex operations, where a complex multiplication and accumulation is considered a single oper- 
ation. As such, we do have not accounted for the 33 percent saving in floating point operations for 
the DFT as described in Section [A. 3. 41 Also, for clarity, we have not included the additional DFTs 
required to support nonlinear dispersion ((5mix) see Section lA.2p . as this substantially complicates 
the analysis, is only significant on the longest baselines and roughly balances the 33 percent saving 
described above. 
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To begin, w e assume the half width beam of a parabohc reflector, at 25% of the peak amplitude 



is (|Cordeslll997l l: 



emax = 0.585-^. (A33) 

The full width of the beam at 25% amplitude is 2^rnax and we set the phase center to the 
center of the primary beam. If we compute only the required DFT bins (as described in Section 
IA.3.ip , the number of DFT bins that must be computed for a single DM over all baselines and the 
full primary beam, is given by: 

M-l M-i ^ 

Nbft = J2 Y1 — 6pg2sin6lmax (A34) 

^ M-l A/-1 

= 2— singmax bpg (A35) 

p=0 ij=0,(j^p 

= 2-sin0max^(M-l)6 (A36) 
c 2 

^ -Oraa^MH^^ (A37) 

c 

where we have employed the small angle formula for sin, and assumed a distribution of baselines 
such that the average baseline length is approximately half of the maximum baseline length. Each 
of these DFT bins requires ~ = -B operations per second (using the block-based sliding DFT. See 
Section IA.3.2P , and assuming we measure A^dm dispersion measures then we require Adm different 
DFT banks. The operations rate for the DFT step is, therefore, given by: 



Cdft = BNbftNbmPjNpoi (A38) 
A B'^ 

^ 0-585;^— ^^'^maxA^DMP/A^pol (A39) 

where Pf = 0.61 is a factor to account for factorizing the DFTs across the DM banks as described 
in Section IA.3.31 and A^poi is the number of polarizations. 

To form an image for a given dispersion measure, a dot product with the truncated amplitude 
response function, across all baselines must be performed for every pixel. The number of pixels in 
an image is: 



iVpix = [^s^j (A40) 
= 1.92 U^^y (A41) 
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For each pixel, we require a dot product with the response function per basehne, implying 



A^ops-per-pixcl = (2F + l)M^ . (A42) 

An image is formed per DM bank at a rate given by T^nt where Tj is the dispersion delay 
associated with the ith DM of interest, and > 1 is the time over sampling factor. If we choose 
set of DM banks that is a geometric progressior0 according to: 



Ti = To{l + ty 0<i<iVDM (A43) 

with e < 1 an overlap factor which can be chosen by a trade-off between computation and SNR. 
The number of DM banks required to cover the ra,ng6 of DJVIs from Tq to ^max 

is given by: 

and, images are produced at a rate: 



log(T^ax/To) 

log(l + e) ^ ^ 



A^DM-l 

= y ^. (A46) 

= ^^DM (A47) 



where 

l-(l + e)-^DM+i 



Finally, the operations rate for the imaging step is: 



(TDM = — / , — (A48) 



Cimg — A^pixA'ops— per— pixclA^imageA'pol 

^ 1.92(K,^?^y(2F + l)M2^aDMArpoi 



^We can choose a geometric progression for the DM bank lengths, and, when Ni > No, round Ni to an integer 
multiple of Nq/L to take maximum advantage from factorization as required in Section TA. 3. 31 
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A. 4- 3. Data rates 

The DFT step takes in voltages for all antennas and produces A^dft outputs sufficient to 
produce images at a rate of A^image- The data rate between the DFT step and the imaging step is 
therefore 

-RoFT-out = ^DFTP/^image-^pol^bytes-pcr-DFT-bin- (A49) 

The data rate at the output of the imaging is: 

^img— out — -^imagc-'%3ix-^pol-^bytes— per— pixel- (A50) 



B. The Chimageator: Analysis and implementation 



In this appendix we describe methods for efficiently implementing The Chimageator. We also 
derive equations for the resolution and operations, and data rate requirements. 



B.l. Non-linear dispersion 

The Chimageator is also affected by non-linear dispersion. As with The Chirpolator (see 
Section IA.2p , the mixing frequency in the non-linear case is no longer constant resulting non-linear 
trajectories in n-k space and higher /co^max- The solutions may be to increase the time oversampling 
Kt and spatial oversampling Kg with additional phase correction to recover coherence. 

As The Chimageator is not really very competitive in the near term (Fig. [6|) we leave a detailed 
treatment for a future paper. 



B.2. Implementation optimizations 

B.2.1. Optimizing operations in the imaging step 

In order to form an image with The Chimageator, we need not sum across all possible trajecto- 
ries. From Fig. [SJ one can see that there is a family of arrival angles and dispersions with the same 
gradient, but whose durations, T, differ. For trajectories on the same gradient, the results for all 
chirp durations can be computed by cumulative sum, i.e. the Px^iOi) can be computed recursively 
from the result -PTi_i(^j) where Tj, Tj-i, 6i and 9j are chosen to have the same gradient in n-k 
space. The gradient of the trajectory is given by: 

fe,(T,^) = l^sin(e) (Bl) 
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The requirement for a shorter trajectory to be calculated from a longer trajectory implies: 



ko{Ti,ei) ~ (B2) 



sm 



^ T, ~ T,^i—^ (B3) 
sm Uj 

where the equivalence of the gradient can be traded depending on SNR and computational require- 
ments. 

Because of this recursive property, there are a much smaller number of independent calculations 
required to search through the DM (equivalent to the Tj) and 6 space than one might naively expect. 
We require only enough operations to calculate the trajectories that end on the rectangle bounded 
by A^max = /sTlnax in the n axis, and ko^max = /s^max sin e^max/c in the k axis, as shown in Figure 
[5j All shorter trajectories can be obtained as partial sums of calculation of the longer trajectory 
with the same (or similar) gradient. 

This optimization works because all trajectories are linear, which means all the shorter trajec- 
tories can be constructed from partial sums of a single longer trajectory with the same gradient . 
When considering the non-linear dispersion, the trajectories are no longer linear and a short tra- 
jectory does not lie along the path of a single long trajectory. One possible solution is to consider 
the trajectories as piece-wise linear. The shorter trajectories can then be constructed from the 
piece-wise partial sums over a number of long trajectories. We leave a detailed treatment of this 
approach to a future paper. 



B. 2. 2. Sampling 

For typical interferometers and dispersions, the gradient of the trajectory is reasonably small, 
which implies an integrate-and-dump operation after the gridding and FFT step can reduce the 
required data volumes and downstream processing requirements. If we assume a trajectory of 
duration T is sampled nt times, and assume the oversampling is proportional to the final DFT bin 
as follows: 

K-t = ^0,end'tt,0 (B4) 

= /,^sin(0)Kt,o (B5) 
then we can produce a sequence of integrated samples indexed by m: 



fsT/Kt 

^fcN = ^ Xk[n + m.fsT/Kt] 

n=0 



(B6) 
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from which we can form the intensity image in Equation [66] in much the same way, but at a reduced 
rate. 

A long integration time will smear out the signal and result in a loss of coherence, which results 
in a practical limit for how small «;t,o can be made. Figure [U] illustrates the amount of coherence 
loss which is achieved for a given value of ki^q. This plot suggests that oversampling the trajectory 
over 5 times during a its duration is sufficient to keep the coherence above 95 per cent. 

B.2.3. Spectral smearing 

As with The Chirpolator, when the instantaneous frequency of the chirp is between discrete 
DFT bins, the energy is spread out along all DFT bins. To recover some SNR in this case, we 
sum along an 2F additional terms in the frequency direction (effectively widening the trajectory) 
to improve the SNR, in much the same way as described in Section 13.41 

B.3. Performance 

B.3.1. Resolution 

The Chirpolator and Chimageator have the same resolution characteristics. This is demon- 
strated by considering ko^max for the two methods. With The Chimageator, k^^max occurs when the 
sample number is the final sample of the chirp, i.e. Ti — f — BT . For The Chirpolator, A^o, maa; 
occurs on the the maximum baseline, hmax- In either case, it has a value: 

A:0,max = S^^^ Sin(^max) (B7) 

and the resolution is given by Equation [A32] 

B.3.2. Operations rates 

The first steps in Chimageator processing are the multiplication and gridding stages. Until 
now, we have assumed a uniform linear array, which makes gridding reasonably straightforward. For 
more complex geometries, a larger gridding kernel is required. A trade between the size of gridding 
support and the quality of the images is outside the scope of this paper, but for dimensioning 
purposes, one can consider a 7 x 7 pixel grid kernel, resulting in A'ops-per-grid-point — 50 including 
the multiplication of the voltages from the two antennas. 

A grid point must be formed from each pair of antennas at the sampling rate, resulting in an 
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operations rate for gridding of: 

M 

C'gridding — fs~^{-^I l)-^ops— per— grid— point (-^^) 

Assuming a spatial oversampling of Kg the number of pixels in the grid plane is: 

r (B9) 

The operations rate for the spatial FFT step is: 

CFFT = /«iVpixlog2iVpix, (BIO) 

and the operations rate for the integration step is: 

Cint = /siVpix. (Bll) 

The total pre-integrator operations rate is therefore: 

C'prc—iiit— total — Cgridding + CpFT + C'int (B12) 

The total pre-integrator operations rate is dominated by the FFT for sparse arrays, while for 
dense arrays, it is dominated by the gridding. 

To compute the data and operations rates of imaging, we begin by assuming we use A'dm 
logarithmically spaced set of trial DMs as described in Equations IA43I and IA441 and that each 
trajectory is sampled at the rate: 

R^ = ^ (B13) 

The computations are then broken into the two types of trajectory shown in Figure [5J the 
trajectories with fixed angle ^max and variable Ti, and the trajectories with fixed time Tmax, and 
variable angle, 9i. In each case, a single trajectory requires K,t{2F + 1) operations per integration 
step. Therefore, the computation rate for a single trajectory is: 



Ctraj = RM2F+1) (B14) 
= ^^^^ (B15) 

= 7^ (B16) 
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We now consider the total operations rate for the one-dimensional case and assuming half the 
beamwidth. 

The fixed angle trajectories have fixed A;o,ond = ^o,max) and variable Tj, resulting in an opera- 
tions rate of: 



Ce^^^ = ^ C'traj (B17) 
i=l 



^o,cnd^t,o(2-^ + 1) (B18) 



=1 * 



= ^0,ma.4o(2i^ + l)E7^ (B19) 

1=1 

= fco'max4o(2^ + l)™ (B20) 

where cjdm is defined in Equation IA48[ 

The fixed time trajectories, have fixed Ti = Tmax and variable A;o,end) resulting in an operations 
rate of operations rate of: 



fc=0 

= 2^ (B22) 

fc=0 * 

fco.max^. (A;/k,)2k2 (2F+1) 

= Yl T 

2 ^0 

= -^^(2F + 1) Y (B24) 

J max 

0,max('^^s ,max + l)(2K,fco,max + l) (B25) 

D-^ max f^s 

^ -4o'^s(2i^ + l)/^o'max (B26) 



ST, 



max 



where Ks is the desired spatial oversampling. The total operations rate is the sum of the two sets 
in the one-dimensional case. In the two dimensional case, the computation is squared, so that the 
total operations rate in 2-dimensions, for the full beamwidth is: 
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C-Totai = (2(C',_ + Ct^))' (B27) 



B.3.3. Data rates 

Assuming this integrate-and-dump operates at the highest rate Rq = Kt/To and the longer 
integrations can be formed from the short integrations in the imaging step, then the data rate at 
the output of the integrate and dump step is: 



R = ATpixA^bytes-per-pix-Ro (B28) 
,max'^s) -kbytes— per— pix^O,max'^t,o/^0 (B29) 
= 4K^Kt,ofeo,max-^bytes-per-pix/Tb (B30) 
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Fig. 2 



, — Taxonomy of approaches to pulsar searching and measurement 
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Fig. 3. — Schematic illustrating the The Chirpolator operating with two antennas. Top left: Two 
linear chirps are received by antennas p and q, with one delayed by r. Bottom left: After taking 
the product of the two voltage time series, the result Xpg has constant frequency over most of the 
duration of the chirp. Bottom right: The DFT of Xpq yields a peak at ko- 
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Fig. 4. — Response of The Chirpolator to a source off the phase center. As the baseUne length 
(bpq) between antennas increases, the position of the peak in the DFT (/cq) increases Unearly, with 
a gradient given by -Bsin(^)/c (see Equation [36|) . The amphtude of the DFT, D]\f{k — ko), is shown 
to illustrate that there is some smearing of the signal around the expected frequency kQ. 
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Fig. 5. — Trajectories of linear chirps with varying durations (T) and angles of arrival (9) after 
gridding and Fourier transforming with The Chimageator. As the sample number (n) increases, 
the peak of the DFT (/cq) increases linearly (see Equation I63p . Two types of trajectory are shown 
with dashed lines: The ^max case, which corresponds a range of dispersion delays, and a single 
arrival angle at the edge of the field of view; and the Tmax case, where each trajectory corresponds 
to the longest dispersion of interest and a range of arrival angles. The thick line is the trajectory 
corresponding to Tjnax/2 and ^max/2, which lies along the T^^x, ^max trajectory and can be therefore 
be computed from the partial sums along th.6 T^max? ^max 

trajectory. 
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Fig. 6. — Data and operations rates for a dispersed pulse survey as a function of algorithm and 
number of antennas. Survey parameters are given in Table [2l Top Panel: data rate between pre- 
integrator and post-integrator steps. Middle panel: operations rate before the integrator. Middle 
panel: operations rate after the integrator. 'Fourier + Tree' and 'Direct + Tree' signify Fourier 
im aging a. n d Dir ect beamforming respectively to form beams, and using tree dedispersion described 
by iTaylorl (jl974l ) for dedispersion. 
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Fig. 7. — Errors in delay calibration reduce the amplitude of the detection metric. Here we plot 
the formation of the detection metric P{9) as the vector (i.e. complex) sum of the phase corrected 
DFT results from three antenna pairs. In the ideal case (red), the phase correction {^pq) perfectly 
corrects for the known phase in the DFT bins {Xpg), and each result has the same absolute phase. 
The resulting detection metric (i-ideai(^))) is fully coherent. If delay errors are present, each DFT 
bin has a residual phase that is different for each antenna pair. The resulting detection metric 
(-forr(^)) has a smaller amplitude, because the vectors do not add into a straight line. 
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Fig. 8. — Fraction of known pulsars that violate the isolated chirp condition, as a function of system 
bandwidth and center frequency (/c). Known pulsars are taken from the ATNF pulsar catalog. 
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Fig. 9. — The simulated mixing frequency as a function of time for a single antenna pair of The 
Chirpolator (see Equation eq:approxerrorl). A range of approximations are shown. The parameters 
for this simulation were: a DM of lOOcm^^pc and a bandwidth of 400 MHz centered at 1.4 GHz, 
9 = 0.5 degree and a baseline of 1 km. 
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Fig. 10. — Error in mixing frequency 6mix as a function of center frequency and bandwidth (B), 
assuming DM of 100cm~^pc, 1 km baseline and 6 = 0.5 degree. 
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Fig. 11. — One method for handling nonhnear dispersion with The Chirpolator is to increase 
the oversamphng rate. Above is the required oversamphng rate (Kt) for a nonhnear chirp vs the 
geometric delay (r) for single a baseline operating at fc = 1.4 GHz, B=400 MHz and a range of 
different coherence losses (rj). The vertical dashed line is the geometric delayfor 9 = 0.5° and a 
baseline of 1 km. The required oversamphng rate is independent of DM. 
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Fig. 12. — The DFT can be efficiently computed in a sliding window manner (see Equation IA17|1 . 
For a DFT bin number k, the current value of the DFT bin (5fc[m]) is formed by taking the previous 
value of the DFT bin (5fc[m — 1]), adding the most recent partial DFT (Vfc[m]) and subtracting the 
oldest partial DFT (Vfc[m - N/L]). 
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Fig. 13. — The calculation of some DFTs can be factorised into the sum of two adjacent DFT 
results. In this example, we illustrate how to calculate the k = 2 bin of the length 2A'^o DFT 
{82'^° [n]) by adding the results of two, adjacent, non-overlapping k = I bins of the length A^o DFT 
{S^° [n]+S^° [n—Nol). In the notation of Section PV.3.31 this example corresponds to k = d = D = 2. 
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Fig. 14. — The Chimageator can operate with non-hnear dispersion as long as a sufficiently high 
oversampling rate is chosen. Above is the required oversampling factor (Kt,o = /^t/^o.end) as a 
function of trajectory gradient {k), to maintain a range of coherence loss levels. The simulated 
array had centre frequency 1.4 GHz, bandwidth 3 MHz, 4 antennas and 500 m spacing. The input 
signal was a linear chirp with dispersion delay corresponding to 20 cm~^ pc. The trajectory gradient 
was calculated for 100 arrival angles from to 3 degrees. The increase at small gradients is due to 
ko,end — for Small gradients but Kt = I. 



